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ABSTRACT 

In  this  paper  we  present  numerical  solutions  to  several  optimal  control  problems  for  an 
unsteady  viscous  flow.  The  main  thrust  of  this  work  is  devoted  to  simulation  and  control 
of  an  unsteady  flow  generated  by  a  circular  cylinder  undergoing  rotary  motion.  By  treating 
the  rotation  rate  as  a  control  variable  we  formulate  two  optimal  control  problems  and  use 
a  central  difference/pseudospectral  transform  method  to  numerically  compute  the  optimal 
control  rates.  Several  types  of  rotations  are  considered  as  potent  ial  controls  and  we  show  that 
a  proper  synchronization  of  forcing  frequency  with  the  natural  vortex  shedding  frequency 
can  greatly  influence  the  flow.  The  results  here  indicate  that  using  moving  boundary  controls 
for  such  systems  may  provide  a  feasible  mechanism  for  flow  control. 
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1.  INTRODUCTION 


The  benefits  received  from  applying  control  mechanisms  in  viscous  flows  has  been  long 
realized  since  the  pioneering  work  of  Prandtl  [30].  Because  of  a  growing  interest  in  controlling 
the  behavior  and  structure  of  fluid  flows,  various  topics  in  flow  control  have  recently  become 
a  subject  of  research  focus  [14,  1,  13,  34,  4,  15,  26,  35,  9,  16,  20,  27,  28,  29,  36],  A  potentially 
important  application  of  flow  control  is  the  enhancement  of  the  aerodynamic  characteristics 
of  future  advanced  aircraft.  Although  recent  research  has  produced  a  variety  of  methods 
to  achieve  desired  performance  and  design  goals  for  maneuvering  vehicles,  it  is  now  widely 
recognized  that  further  gains  will  most  likely  come  from  the  application  of  various  types 
of  flow  control  mechanisms  [11,  12].  Consequently,  the  utility  of  flow  control  becomes  a 
critical  issue  in  the  design  process  which  may  provide  real-time  effect  for  many  important 
applications,  such  as  highly  instantaneous  maneuvers  for  the  super-maneuverable  aircraft 
[18],  and  the  optimum  design  of  aerodynamic  configurations  [21].  Considerable  effort  has 
been  devoted  to  the  improvement  of  control  mechanisms.  HoweveT,  the  principal  progress  to- 
date  has  been  essentially  accomplished  by  experimental  investigations,  w'liile  most  analytical 
and  numerical  approaches  have  remained  in  its  infancy  due  to  the  complexity  of  the  problems. 

One  of  the  most  practical  applications  of  control  mechanisms  in  flow  systems  has  been 
boundary-layer  separation  control.  In  this  area  of  research,  several  methods  have  been  de¬ 
veloped  experimentally  to  provide  various  control  mechanisms,  for  example  moving  surfaces, 
blowing,  suction,  injection  of  a  different  gas,  etc  [31,  32].  In  particular,  it  has  been  demon¬ 
strated  in  a  number  of  experiments  by  Modi  et  al.  [25,  24]  that  moving  surfaces  can  effectively 
provide  boundary-layer  control.  In  their  experiments,  the  boundary-layer  flow  is  controlled 
by  an  application  of  two  rotating  cylinders  located  at  the  leading  and  trailing  edges  of  an 
airfoil.  It  has  been  shown  that  this  mechanism  can  prevent  flow  separation  by  retarding  the 
initial  growth  of  the  boundary  layer,  with  the  important  consequences  of  lift  enhancement 
and  stall  delay.  For  instance,  when  the  speed  ratios  (which  represents  the  ratio  of  cylinder 
surface  speed  to  the  freestream  speed)  of  both  cylinders  were  set  at  a  constant  value  of  4,  the 


1 


results  indicated  a  200%  increase  of  the  maximum  lift  coefficient  compared  with  the  reference 
airfoil  (in  which  no  rotating  cylinder  is  attached).  In  spite  of  the  fact  that  considerable  aero¬ 
dynamic  benefits  were  gained  by  changing  the  cylinder  speed  ratio,  in  their  experiments  the 
speed  of  rotation  was  performed  merely  with  constant  values.  However,  it  should  be  noted 
that  if  the  rotating  cylinder  mechanism  is  applied  to  a  region  of  unsteady  flow,  a  constant 
rotation  rate  may  not  correspond  to  the  optimal  performance  when  an  airfoil  is  undergoing  a 
rapid  maneuver.  This  type  of  result  provided  the  motivation  for  us  to  consider  a  fundamen¬ 
tal  problem  regarding  unsteady  flow  control  by  means  of  a  time-dependent  moving  surface 
mechanism.  In  order  to  keep  the  problem  reasonable  and  yet  practical,  we  selected  a  model 
for  the  numerical  study  of  controlling  the  temporal  development  of  the  flow  field  around  a 
rotating  cylinder. 

The  most  distinguishing  feature  of  a  rotating  body  traveling  through  a  fluid  is  that  it 
experiences  a  transverse  force  acting  in  a  direction  perpendicular  to  that  of  flowing  stream 
[37],  In  the  past  few  decades,  research  on  the  problem  of  a  uniform  stream  past  a  cylindrical 
rotating  body  has  been  the  subject  of  many  experimental  and  numerical  investigations.  See 
the  papers  by  Taneda  [39],  Mo  [23]  and  Tokumaru  and  Dimotakis  [40]  for  a  cylinder  under¬ 
going  rotary  oscillations,  Prandtl  [31],  Taneda  [38],  Koromilas  and  Telionis  [22],  Coutanceau 
and  Menard  [8],  Badr  and  Dennis  [3],  Badr  et  al.  [2]  and  Chen  et  al.  [7]  for  a  cylinder  with 
a  constant  speed  of  rotation.  However,  most  of  these  results  are  primarily  focused  on  the 
study  of  formation  and  development  of  vortices  in  a  cylinder  wake.  It  appears  that  the  effect 
of  the  rotation  rate  on  the  cylinder  forces  exerted  by  the  fluid  has  received  far  less  attention, 
despite  the  fact  that  it  has  many  important  practical  engineering  applications. 

The  main  thrust  of  the  current  investigation  is  on  simulation  and  control  of  an  unsteady 
flow  generated  by  a  circular  cylinder  undergoing  a  combined  (steady  or  unsteady)  rotary 
and  rectilinear  motion.  By  treating  the  rotation  rate  as  a  control  variable  in  this  model, 
we  consider  several  problems  concerning  the  temporal  development  of  forces  on  a  rotating 
cylinder  in  response  to  a  variety  of  time-dependent  rotation  rates.  The  computational  results 
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provide  considerable  insight  into  the  problem  of  controlling  such  flows. 

2.  MATHEMATICAL  AND  NUMERICAL  FORMULATIONS 

In  this  section  the  governing  equations  and  the  particular  numerical  method  used  in  this 
work  are  described.  We  consider  control  problems  for  a  two-dimensional  viscous  incompress¬ 
ible  flow  generated  by  an  impulsively  started  circular  cylinder.  The  cylinder  is  translated 
with  a  constant  rectilinear  speed  U  normal  to  its  generator  and  is  simultaneously  rotated 
with  a  time-dependent  angular  velocity  Ct(t)  about  its  axis.  Although  there  are  various 
formulations  and  numerical  techniques  for  the  solving  of  steady  and  unsteady  flow  past  a 
rotating  cylinder  [19,  3,  23],  in  this  work  the  problem  is  investigated  numerically  by  solving 
a  velocity/vorticity  formulation  of  the  Navier-Stokes  equations  with  an  implementation  of 
the  Biot-Savart  law.  The  numerical  approach  used  in  the  present  study  is  the  one  developed 
by  Chen  [5]  for  the  problem  of  a  circular  cylinder  oscillating  in  a  rectangular  box.  It  is  based 
on  an  explicit  finite-difference/pseudo-spectral  technique  to  yield  time  accurate  solutions 
to  the  governing  equations.  This  numerical  algorithm  was  further  modified  to  investigate 
an  unsteady  flow  around  a  rotating  cylinder  undergoing  various  constant  rotational  speeds 
[7,  28],  and  time-dependent  rotation  rates  [26,  27]. 

Several  prominent  features  in  the  application  of  an  integral  representation  for  flow  kine¬ 
matics  are  emphasized  in  §2.2.  This  integral  method  proposed  by  Wu  and  Thompson  [46] 
provides  the  basic  link  between  the  velocity  and  vorticity  fields  throughout  the  numerical 
procedure.  The  boundary  vorticity  at  the  solid  surface  can  be  easily  calculated  by  the  ap¬ 
plication  of  this  integral.  Moreover,  unlike  other  numerical  approaches,  the  imposition  of 
the  artificial  far-field  boundary  condition  for  the  velocity  is  not  necessary  once  the  vorticity 
values  are  known  everywhere  in  the  domain  of  interest. 

2.1  Governing  Equations 

In  the  velocity/vorticity  formulation  of  the  Navier-Stokes  equations,  the  governing  equa¬ 
tions  consist  of  the  vorticity  transport  equation  and  the  vector  Poisson  equation  for  the  veloc- 
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ity.  Thus,  for  a  two-dimensional  unsteady  viscous  flow  in  incompressible  fluid,  the  Cartesian 
coordinate  form  of  the  governing  equations  for  the  vorticity  and  velocity  can  be  expressed 
in  the  dimensionless  form  as 

%  +  (1) 

and 

V2u  =  -Vx(  we,),  (2) 

where  u  is  the  velocity  field,  u>  is  the  vorticity  field,  and  ez  is  the  unit  vector  in  the  direction 
of  ^-direction.  All  the  variables  are  made  dimensionless  by  means  of  the  characteristic 
quantities.  The  cylinder  radius  a  is  used  as  the  length  scale  while  ajU  is  used  as  the  time 
scale.  The  Reynolds  number  Re  =  2 Uajv  is  based  on  the  cylinder  diameter  2 a  and  the 
magnitude  U  of  the  rectilinear  velocity. 

A  non-rotating  reference  frame,  translating  with  the  cylinder  is  employed.  In  this  frame 
the  dimensionless  boundary  conditions  for  the  problem  of  a  rotating  cylinder  (with  a  time- 
dependent  angular  velocity  O(t)e^)  can  be  written  as 

u  =  -a(t)yex  +  a(t)xey  for  (x,  y)  €  T  (3) 

and 

u  —  ex  for  yjx2  -f  y2  — >  oo,  (4) 

where  F  denotes  the  impermeable  solid  boundary  of  the  cylinder.  The  angular/rectilinear 
speed  ratio  a(t)  =  U(t)a/U  is  the  primary  control  parameter  throughout  this  work. 

In  many  practical  numerical  simulations  for  the  laminar  motion  of  viscous  incompressible 
fluid,  both  the  exterior  and  interior  flow  problems,  the  formulation  based  on  velocity/vorticity 
variables  would  provide  some  advantages  over  the  primitive-variable  formulation.  This  ve¬ 
locity/vorticity  formulation  is  especially  well  suited  to  treating  initial  development  of  the 
flow  generated  by  impulsively  started  bodies,  in  which  a  relatively  small  vortical  viscous  re¬ 
gion  is  embedded  in  a  much  larger  inviscid  potential  flow.  Consequently,  the  computational 
domain  may  be  restricted  to  a  smaller  region  where  all  vorticity  contributions  are  contained. 


4 


Furthermore,  the  decoupling  of  the  overall  problem  into  its  kinematic  and  kinetic  aspects  of 
the  flow  field  gives  an  additional  convenience.  Throughout  this  work  the  vorticity  transport 
equation  described  above  may  be  viewed  as  kinetic  process  of  the  flow  field  in  which  the 
distribution  of  vorticity  is  interplayed  by  the  process  of  convection,  as  well  as  the  effects  of 
viscous  diffusion. 

2.2  Integral  Representation  for  Flow  Kinematics 

Having  defined  the  governing  equations  for  the  rotating  cylinder  problem  described 
above,  we  can  now  examine  several  advantages  of  an  integral  formulation  used  in  this  work. 
In  the  formulation  of  exterior  flow  problems,  it  is  well  known  that  one  of  the  main  numerical 
difficulties  is  related  to  the  proper  imposition  of  a  prescribed  limit  value  at  the  far-field  for 
the  unbounded  physical  domain  in  which  the  flow  takes  place.  It  should  be  pointed  out  that 
in  any  numerical  simulation  one  has  to  restrict  the  exterior  infinite  domain  to  be  finite  with 
an  artificial  boundary.  However,  the  far-field  boundary  in  (4)  fails  to  represent  the  exact 
characteristics  at  the  outer  perimeter  of  the  finite  computational  domain.  Therefore,  in  many 
practical  applications,  instead  of  directly  applying  the  far-field  boundary  condition  (4)  one 
often  tries  to  avoid  the  difficulty  by  utilizing  various  asymptotic  boundary  conditions  at  large 
distance  (e.g.  the  application  of  potential  flow  or  Ossen  expansion),  while  others  introduce 
a  mapping  of  the  infinite  domain  onto  the  finite  one  by  means  of  a  suitable  coordinate 
transformation. 

The  other  difficulty  encountered  in  the  simulation  of  viscous  flow  is  that  of  prescrib¬ 
ing  the  appropriate  non-velocity  boundary  conditions  at  the  solid  surface.  In  general,  the 
prescribed  pressure  at  the  solid  boundary  is  needed  in  the  application  of  primitive-variable 
(pressure/velocity)  formulation,  while  the  boundary  vorticity  is  required  for  the  formulation 
based  on  the  velocity/vorticity  (or  stream-function/vorticity)  variables.  In  order  to  over¬ 
come  these  two  difficulties,  we  pose  the  kinematic  relation  on  the  problem  by  introducing  a 
general  Biot-Savart  induced  law  described  below. 


In  the  problem  with  a  viscous  fluid,  if  the  velocity  distribution  of  u  are  given,  then  the 
vorticity  field  3  could  be  evaluated  through  the  kinematic  relation  between  u  and  3  described 
by  3  =  V  x  u  and  V  -  u  =  0.  On  the  other  hand,  the  vector  Poisson  equation 

V2u  =  —V  x  3,  (5) 


obtained  from  the  continuity  equation  and  the  definition  of  vorticity,  can  be  used  to  determine 
the  velocity  field  from  a  given  vorticity  field. 

For  a  general  viscous  flow  in  a  region  D  bounded  by  an  inner  boundary  F  and  outer 
boundary  F',  we  can  recast  the  kinematic  part  of  the  problem  into  an  equivalent  integral 
formulation 


u(ro,<) 


-i  /  /  rl IdA 

c  J  Jd  |r  —  rop 

_  1  f  K(r,Q  •  n](r-  fQ)  -  [ ub(r,t )  x  n]  x  (r-  r0) 
c  Jr+r<  |r  —  r0|d 


(6) 


where 

47r,  for  d  =  3 

2rr,  ford  =  2.  ('j 

In  the  above  integral  representation,  tq,  is  the  boundary  velocity,  n  is  the  outward  normal 
unit  vector  and  d  is  the  spatial  dimension.  The  subscript  u0”  denotes  the  field  point  where 
the  velocity  field  is  evaluated. 

By  applying  the  no-penetration  and  no-slip  conditions  to  the  rotating  cylinder  problem 
considered  in  this  study,  in  two  dimensions  equation  (6)  can  be  written  in  terms  of  the 
rectilinear  velocity  Uex  and  the  angular  velocity  of  the  solid  body  D  which  is  known 

as  the  generalized  Biot-Savart  law  of  induced  velocity: 

t)  x  (f—  ro)  ,  . 


For  detailed  discussions  of  the  integral  representation  for  viscous  flows  the  reader  is  referred 
to  [46,  45].  Equation  (8)  represents  the  kinematic  relationship  between  the  velocity  and  vor¬ 
ticity  fields  of  an  infinite  domain  which  is  expressed  as  an  integral  form.  The  first  integral 


|f-ro|2 
*,t)  x  (r-rp) 
|r  —  ro|2 


dA  +  U. 


(8) 
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represents  the  contribution  of  the  vorticity  field  to  the  development  of  the  velocity  field  over 
the  region  D  occupied  by  the  viscous  fluid,  while  the  second  integral  gives  the  contribu¬ 
tion  from  the  rotation  of  th'1  ^ud  body.  Notice  that  both  the  no-penetration,  the  no-slip 
boundary  conditions  on  the  solid  surface  and  the  far-field  boundary  conditions  are  implicitly 
imposed  within  the  integral  (8).  To  be  precise,  if  we  let  the  calculated  field  point  r0  go  to 
infinite  in  (8),  the  vorticity  contribution  v’ill  approach  zero.  Consequently,  the  boundary 
condition  at  infinity  is  satisfied  exactly  in  this  integral  representation  of  the  kinematic  rela¬ 
tion.  This  indicates  that  the  difficulty  resulting  from  the  imposition  of  far-field  condition  is 
removed  by  an  application  of  (8). 

In  numerical  simulations,  it  has  been  reported  that  the  imposed  boundary  conditions  at  a 
large  distance  from  the  body  will  significantly  influence  the  accuracy  of  the  overall  numerical 
solution  [10,  19].  It  is  therefore  necessary  to  devise  a  technique  to  impose  the  appropriate 
condition  at  the  outer  boundary  of  the  computational  domain.  The  integral  approach  in 
(8)  provides  a  useful  method  to  accomplish  this.  To  be  precise,  the  integral  representation 
(8)  permits  us  to  determine  the  velocity  (point-by-point)  explicitly  if  all  vorticity  values  are 
known  everywhere  in  the  domain  of  interest.  It  will  exhibit  a  more  realistic  behavior  at  the 
outer  perimeter  of  computational  domain  than  those  asymptotic  techniques  employed  by 
other  formulations.  Namely,  if  the  computational  domain  is  large  enough  to  contain  all  of 
vorticity  generated  around  the  solid  boundary  prior  to  a  certain  time,  then  at  this  instant 
the  velocity  on  the  outer  perimeter  of  computational  domain  can  be  evaluated  directly  by 
the  numerical  integration  of  (8)  with  all  the  known  vorticity  field  in  the  domain. 

Additionally,  if  we  apply  equation  (8)  to  the  points  of  r*,  on  the  solid  boundary,  then  the 
integral  formula  becomes 


v(r,t)  x  (f-  rb) 
|r  -  rl|2 

2Q(r,£)  x  (r  -  r^) 
jr  -  rb |2 


dA 

dA  +  U. 


(9) 


The  boundary  vorticity  values  are  contained  in  this  integral.  Hence,  by  using  the  prescribed 
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motions  of  body  (i.e.  U  and  { l(t))  and  all  of  the  known  vortieity  field  in  the  domain,  the 
determination  of  vortioity  values  at  grid  points  on  the  solid  boundary  can  be  achieved  by  the 
numerical  integration  of  equation  (9).  However,  many  implementations  of  flow  simulations 
make  use  of  approximate  formulae  for  the  boundary  vortieity,  which  may  lead  to  excessive 
error  and  thereby  destroying  the  accuracy  of  the  solutions. 

When  we  solve  the  velocity  field  in  an  infinite  domain  we  need  only  to  take  the  viscous 
region  (J  ^  0)  into  account,  while  contributions  from  the  inviscid  region  vanishes  (T  =  ()}. 
Hence,  an  application  of  this  integral  representation  allows  one  to  confine  the  computation 
to  the  viscous  region  corresponding  to  the  non-negligible  vortieity  portion  of  the  flow,  in 
contrast  to  employing  domain  methods  (i.e.  the  finite-difference  and  finite-element  met  hods) 
in  which  both  potential  and  viscous  regions  are  needed  in  a  simultaneous  solution.  This 
implies  the  computational  domain  can  be  significantly  reduced  and  lienee  a  smaller  number 
of  grid  points  are  needed  than  those  required  by  standard  finite-difference  and  finite-element 
methods. 

2.3  General  Body-Fitted  Coordinate 

In  order  to  accommodate  problems  with  a  time-dependent  domain  in  the  physical  (r,  y) 
space,  the  vortieity  transport  equation  (1)  and  Poisson  equation  (2)  can  be  recast  in  terms 
of  the  time- varying  generalized  body-fitted  coordinate  system  (£,  r/)  as 

xt , 

—  -J\^y  fj  — 

-J  [yr,(UUj)i  -  +  Xn(VU,)t  -  XdVUJ)r>] 

2  2 

+  ~  +  7^)  +  —  ( (10} 

and 

oua  -  2fhi0l  4-  7 u,jn  -f  P(PH  +  Quv)  =  .7(xr,u.v  -  x^<v), 

2/P  {tj  4  7  T  J  (Pv^  4  QUjj)  —  J (f/rj-'-V 

where 

f  o  ~  x'*  +  y*,  p  =  Tf-r^  4-  yf.yn,  7  =  x2(  4-  t/|, 

\  P  T  £yy,  Q  —  7jr],  1]yy  ,  ./  =  J" ,,  ~  ■ 


y^)  ~  jiUtXq-UqXt) 


(12) 


Here  J  is  the  Jacobian  of  the  coordinate  mapping  between  the  physical  space  (,r,y)  and  the 
computational  space  (f,  jj).  All  subscripts  in  ( 10)-(  12)  denote  partial  differentiation.  Also, 
u  and  v  denote  as  the  velocity  component  in  x-  and  {/-direction,  respectively.  Notice  that 
this  coordinate  transformation  will  introduce  two  additional  time-dependent  terms  rt  and 
yt  when  a  time-dependent  flow  domain  in  the  physical  space  is  considered.  However,  one 
can  calculate  the  solution  in  a  fixed  time-independent  computational  grid  in  this  generalized 
coordinate  system  because  of  the  inclusion  of  xt  and  y,  in  the  governing  equations.  Therefore, 
even  when  a  time-dependent  domain  is  considered,  the  interpolation  of  boundary  conditions 
and  grid  points  employed  in  the  physical  space  formulation  is  not  necessary  in  this  generalized 
coordinate  system.  Moreover,  it  is  also  convenient  to  distribute  the  grid  points  within  the 
domain  of  interest  if  a  complex  flow  geometry  is  encountered.  In  fact,  such  a  distribution  of 
grid  points  is  particularly  important  for  simulations  in  viscous  flows,  where  the  grid  points 
clustering  near  the  surface  is  necessary  to  resolve  the  large  gradients  which  appeared  in  a 
thin  region  due  to  the  viscous  eflfec* 

Decause  we  wish  to  extend  the  results  presented  here  to  more  general  problems  with 
time-dependent  domains,  the  computer  code  was  written  under  this  general  body-fitted 
coordinate  formulation.  However,  in  this  work,  the  polar  coordinate  (r,  9)  is  used  for  the 
study  of  an  unsteady  flow  around  a  rotating  cylinder  and  the  flow  domain  considered  here 
is  time-independent  since  a  reference  frame  translating  with  the  cylinder  is  used.  If  a  body 
of  arbitrary  shape  other  than  circular  cylinder  is  considered  (e.g.  a  rotating  elliptic  cylinder 
or  airfoiD  then  the  domain  of  interest  becomes  a  time-dependent  feature  even  though  the 
reference  frame  is  translating  with  the  body.  This  occurs  because  of  the  rotational  motion 
of  an  asymmetric  bodies  with  respect  to  a  non-rotating  reference  frame. 

2.4  Computational  Procedure 

The  computational  domain  is  discretized  in  the  body-fitted  coordinate  system  (£.  tj)  In¬ 
setting  &  =  iAC  t) j  =  jArj,  and  i  -  1 ,  -  - ,  A/ ,  j  =  1,  — ,  Ar,  where  M  and  N  denote  the 
number  of  grid  points  in  the  £-  and  {/-direction,  respectively.  The  numerical  approximation  of 


9 


the  vorticity  at  the  grid  point  (Zi,T]j)  and  time  tn  —  nAt  is  denoted  by  —  un{iA£,jAtj). 
The  vorticity  transport  equation  (10)  is  first  discretized  by  a  second  order  central  difference  in 
the  radial  direction  and  a  pseudospectral  transform  method  in  the  circumferential  direction 
for  all  spatial  derivatives.  This  semi-discretization  form  of  equation  (10),  consisting  of  a 
system  of  ordinary  differential  equations  in  time  can  be  written  as 

—  F(u),  U  —  (o>2,2>  '  ’  •  »k>A/-l,/V-l)T,  (13) 


for  all  the  interior  grid  points.  Therefore,  the  calculation  procedure  to  advance  the  solution 
for  any  given  time  increment  can  be  summarized  as  follows: 

Step  1:  Internal  vorticity  over  the  fluid  region  at  each  interior  field  point  is  calculated  by 
solving  the  discretized  vorticity  transport  equation.  An  explicit  second-order  rational  Runge- 
Kutta  marching  scheme  based  on  the  work  of  [4 1]  is  used  to  advance  in  time  for  (IS) 

The  discretization  of  (13)  in  time  thus  can  be  written  as 

2$*i  ($*1,5*3)  “  $3($i, 9i) 


U’n  +  1  =  <5”  + 


($*3,5*3) 


(14) 


with 

'  $'1  =  F(u>n)At 

4  g~2  -  F(un  +cgi)At  (15) 

.  $3  =  (1  ~  6)  5*1  -  bg2 

where  (:, :)  denotes  the  scalar  product.  In  order  to  ensure  the  stability  of  the  above  nonlinear 
explicit  scheme,  the  two  constants  b  and  c  in  (15)  must  satisfy  be  =  —0.5.  In  particular, 
b  =  —  1  and  c  =  0.5  are  used  in  our  computations.  Although  this  method  is  explicit  in 
nature,  it  may  become  unconditionally  stable  by  the  suitable  choice  of  the  constants  in  (15) 
[17].  In  addition,  this  particular  scheme  allows  one  to  use  a  larger  time-step  than  that  of  the 
three-step  Adams-Bashforth  scheme  used  by  Chen  [5].  This  step  consists  of  the  kinetic  part 
of  the  computational  loop. 

Step  2:  Using  known  internal  vorticity  values  at  all  the  interior  grid  points  from  step 
1,  the  generalized  Biot-Savart  law  of  induced  velocity  (9)  is  used  to  update  the  boundary 
vorticity  values  at  all  the  surface  nodes. 
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Step  3:  At  this  stage,  all  the  vorticity  values  in  the  computational  domain  are  known  at 
the  new  time  level.  Then,  the  velocity  at  points  on  the  outer  perimeter  of  the  computational 
domain  is  calculated  by  the  integral  kinematic  constraint  (8). 

Notice  that  at  each  time  step,  the  numerical  integration  of  the  first  integral  over  the  fluid 
region  in  (8)  is  carried  out  by  means  of  an  isoparametric  formulation  which  is  used  extensively 
in  the  finite  element  method,  while  the  second  integral  can  be  evaluated  analytically  over 
the  solid  body  B.  The  vorticity  distribution  over  each  distorted  quadrilateral  element  in 
the  physical  space  are  actually  performed  over  a  square  in  the  isoparametric  space.  Further 
details  of  the  integration  method  can  be  found  in  Chen  [5]. 

It  is  worthy  of  note  that  the  evaluation  of  the  integral  (8)  involves  a  problem  associated 
with  the  non-uniqueness  of  the  solution.  The  principle  of  vorticity  conservation  imposed 
by  Wu  [45]  resolve  such  difficulty.  For  flows  past  single  or  multiple  solid  bodies,  a  more 
immediate  improvement  to  the  principle  of  vorticity  conservation  is  provided  by  Chen  [5]. 

Step  4-'  The  new  internal  velocity  field  can  be  established  by  solving  the  Poisson  equations 
(11)  with  prescribed  solid  boundary  conditions  and  outer  boundary  conditions  of  the  velocity 
that  have  been  determined  from  step  3. 

The  final  form  of  the  discretized  Poisson  equations  can  be  written 

{ir=l  <'6) 

where  u  and  v  are  two  vectors  of  unknown  interior  nodal  values.  Also,  A  is  a  11-bandcd 
matrix,  while  fj  and  f 2  are  vectors  associated  with  the  known  forcing  terms  and  boundary- 
conditions.  The  resulting  11 -banded  matrix  equations  are  then  solved  by  a  preconditioned 
biconjugate  gradient  routine  [6],  This  step  completes  the  computational  loop  for  each  time 
level. 

One  further  important  point  to  be  noted  in  the  integral  approach  is  that  the  initial  flow 
field  can  be  determined  by  the  same  solution  procedure  described  above  (from  step  2  to 
step  4).  To  be  more  precise,  at  time  t  =  0+  the  unknown  boundary  vorticity  values  is 
determined  by  solving  integral  (9)  with  zero  vorticity  values  everywhere  away  from  the  solid 
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boundary  surface.  This  is  based  on  the  fact  that  the  vorticity  will  concentrate  on  the  body 
surface  in  the  form  of  a  vortex  sheet  immediately  after  the  body  impulsively  started.  This 
approach  therefore  reduces  equation  (9)  to  a  boundary  integral,  indicate  that  the  method 
can  be  viewed  as  the  boundary  element  method  utilized  in  the  problem  of  the  potential  flow. 
Once  the  boundary  vorticity  values  are  obtained,  the  initial  velocity  field  can  be  determined 
by  solving  equation  (11)  with  the  known  velocity  values  at  all  points  on  the  outer  perimeter 
of  the  computational  domain  which  have  been  calculated  by  the  integral  (8).  In  contrast  to 
the  special  technique  used  by  other  methods,  this  integral  approach  enables  the  numerical 
code  to  generate  the  initial  velocity  field  simply  by  the  implementation  of  one  cycle  of  a 
solution  procedure  rather  than  employing  any  additional  treatments. 

3.  RESULTS  AND  DISCUSSIONS 

In  this  section  we  apply  the  numerical  scheme  described  above  to  simulate  and  control 
the  unsteady  flow  around  a  rotating  cylinder  that  undergoes  a  variety  of  steady  and  unsteady 
angular/rectilinear  speed  ratios  at  a  Reynolds  number  of  200.  The  choice  of  this  particular 
Reynolds  number  is  not  due  to  the  limitation  of  the  numerical  algorithm,  it  is  mainly  for 
the  purpose  of  comparing  with  the  existing  experiments  of  Coutanceau  and  Menard  in  the 
case  of  constant  rotation  for  a  rotating  cylinder  [8].  In  this  model,  the  rectilinear  velocity  is 
fixed  as  a  constant  value  while  the  angular  velocity  is  treated  as  a  control  variable. 

Although  the  choice  of  time-dependent  rotation  rates  that  may  be  used  to  control  the 
rotating  cylinder  are  unlimited,  the  computational  results  presented  here  are  restricted  to 
the  following  three  types  of  rotation: 

1.  Constant  speed  of  rotation:  a(t )  =  constant. 

2.  Time-harmonic  rotary  oscillation:  a(f)  =  AsinirFt. 

3.  Time-periodic  rotation:  a(t)  =  A\  sin  ir(F/2)t\. 

Again,  this  choice  was  made  because  it  allowed  us  to  compare  experimental  and  numerical 
results  and  it  matched  the  control  experiments  of  Modi  [25,  24].  All  variables  are  normalized 
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to  the  nondimensional  forms  iu  the  formulation.  In  a  non-rotating  frame  attached  to  the 
cylinder,  the  configurations  for  the  different  controls  considered  in  the  physical  space  are 
sketched  in  Figure  1,  together  with  the  corresponding  time  evolution  of  the  angular  velocity. 

Concerning  the  use  of  a  time-harmonic  rotary  oscillation,  it  is  common  to  define  the 
motion  of  angular  rotation  0(f)  as 

0(f)  =  -dcos2irft,  (17) 

where  6  is  the  angular  amplitude  and  /  is  the  forcing  frequency  of  the  oscillation.  Thus,  the 
associated  time-dependent  speed  ratio  is  given  by 

a(t)  =  —  As'mnFt,  (18) 

where  F  =  2af  /U  is  the  reduced  forcing  frequency  and  A  =  nF9  is  the  normalized  maximum 
rotation  rate  of  the  forcing  oscillation.  In  order  to  attain  high  lift  and  reduced  drag,  previous 
work  for  constant  rotation  rates  [7]  lead  us  to  consider  one  special  type  of  time-periodic 
rotation.  That  is,  the  cylinder  under  control  is  rotated  in  the  counterclockwise  direction 
about  its  axis  with  a  time-periodic  speed  ratio  given  by 

a(t)  =  A\s\mt(F /2)t\.  (19) 

Here,  the  reduced  forcing  frequency  of  this  particular  time-periodic  rotation  is  F.  This  partic¬ 
ular  type  of  rotation  is  expected  to  provide  a  substantial  lift  enhancement  and  drag  reduct  ion 
through  a  proper  choice  of  both  the  angular  amplitude  (thus  the  normalized  maximum  ro¬ 
tation  rate  A)  and  forcing  frequency  (thus  the  reduced  frequency  F).  This  improvement  can 
be  demonstrated  by  comparing  its  respective  force  performances  against  the  time-harmonic 
rotary  oscillation. 

The  major  goal  of  this  paper  is  to  study  the  effect  of  rotation  rate  control  upon  the 
lift  and  drag  on  the  cylinder  surface.  Hence,  in  the  following  discussion  we  concentrate  on 
various  issues  concerning  the  development  of  temporal  forces.  In  a  viscous  flow,  it  is  well 
known  that  the  total  lift  and  drag  forces  are  contributed  by  the  pressure  and  skin  friction  due 


13 


to  the  viscous  effects.  An  important  consequence  of  using  the  velocity /vorticity  formulation 
is  that  the  forces  can  be  directly  evaluated  from  the  known  vorticity  on  the  cylinder  surface. 
Hence,  for  known  vorticity  values  on  the  cylinder  surface,  the  lift  and  drag  coefficients  can 
be  calculated  in  the  r-d  coordinates  by 

Cc(t)  =  CLp(t)+Cif(t)  =  cos  6d9  +  ~  u>(t)r  cos  Odd,  (20) 

and 

CD(t)=CDp(t)+CDf(t)  =  ~£r  sin  OdB  ~  —  j**  u(t)r  sin  ddQ,  (21) 

where  the  subscript  T  denotes  quantities  evaluated  on  the  cylinder  surface.  The  subscripts  p 
and  /  represent  the  contribution  from  pressure  and  skin  friction,  respectively.  In  particular, 
we  denote  the  positive  values  of  Ci  in  the  negative  {/-direction  (as  noted  in  Figure  1). 

To  assess  the  accuracy  of  the  numerical  algorithm,  computations  were  first  performed 
over  a  wide  range  of  constant  speed  ratios  up  to  3.25  at  a  Reynolds  number  of  200.  Several 
particular  speed  ratio  parameters  were  chosen  to  allow  for  the  comparison  against  the  ex¬ 
perimental  work  of  Coutanceau  and  Menard  [8].  Speed  ratios  greater  than  2  are  important 
in  the  study  of  the  possibility  of  suppressing  vortex  shedding  by  an  application  of  higher 
rotation  rates.  The  details  of  the  work  using  constant  rotation  rates  were  reported  in  [7], 
where  several  numerical  solutions  were  compared  and  demonstrated  to  be  in  good  agreement 
with  experimental  results. 

In  the  computations  below,  a  fixed  flow  domain  is  used  and  its  extent  is  essentially 
determined  by  the  time-span  under  investigation.  Namely,  a  larger  computational  domain 
will  be  needed  for  a  longer  time  of  observation.  Here,  a  circle  of  radius  r  =  24  and  r  =  36 
are  chosen  for  the  time-span  of  0  <  t  <  24  and  0  <  t  <  36,  respectively.  These  long  time 
histories  were  necessary  to  demonstrate  the  periodicity  of  the  flow  pattern  and  evolution  of 
the  forces.  When  a  larger  domain  is  considered,  the  mesh  is  increased  in  the  radial  direction 
in  order  to  properly  discretize  the  domain.  A  uniformly  spaced  M  =  128  in  circumferential 
direction  is  used  for  all  computations,  while  N  =  120  and  N  =  180  stretched  grid  lines 
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in  the  radial  direction  are  used  for  the  time-span  under  consideration.  In  the  numerical 
calculations,  small  initial  time  steps  are  taken  in  order  to  contain  the  inherent  numerical 
error  caused  by  resolving  of  the  infinitesimal  vorticity  layer  at  t  =  0+.  A  detailed  discussion 
of  the  algorithm  and  the  accuracy  of  the  initial  flow  field  can  be  found  in  [7], 

3.1  Constant  Speed  of  Rotation 

Figures  2(a,b)  show  calculated  instantaneous  streamline  plots  for  a  constant  value  of 
speed  ratio  of  a  =  2.07  at  Re  =  200  and  compare  these  plots  with  the  experimental  work 
of  Coutanceau  and  Menard  [8].  In  the  computation,  the  non-rotating  reference  frame  is 
translating  with  the  cylinder  while  the  camera  in  the  experiment  is  moving  with  the  cylinder 
as  well.  Excellent  agreement  is  obtained,  despite  the  fact  that  a  high  velocity  gradient  is 
induced  in  the  near  wake  due  to  the  cylinder  rotation.  Although  not  shown  here,  one  partic¬ 
ular  interesting  feature  is  the  difference  between  the  experimental  work  and  our  calculated 
observation  regarding  the  conclusion  of  suppressing  of  vortex  shedding  at  high  speed  ratios. 
In  the  computation  using  high  speed  ratios  (at  Re  =  200,  as  reported  in  [7]),  the  calculated 
equi- vorticity  contours  seem  to  imply  that  vortex  shedding  continues  to  occur  even  at  high 
rotation  rates  (c*  >  2.07).  However,  at  these  high  a,  the  observed  formation  of  the  vortex 
street  behind  a  rotating  cylinder  seems  to  contradict  the  experimental  conclusion  described 
in  [8].  This  difference  is  due  to  the  fact  that  the  experimental  apparatus  was  such  that  only 
10  dimensionless  time  units  of  data  could  be  collected  and  in  part  by  the  flow  visualization 
techniques  used  in  their  experiments.  On  the  other  hand,  it  is  important  to  note  a  recent 
investigation  by  Badr  et  al.  [2]  regarding  the  issue  of  suppressing  of  vortex  shedding.  Their 
tests  were  performed  both  experimentally  and  numerically  at  Reynolds  numbers  of  Re  =  103 
and  Re  =  104.  For  a  rotation  rate  of  a  =  3  at  Re  =  103,  they  show  that  no  other  eddy  is 
created  after  the  shedding  of  two  vortices.  In  addition,  the  temporal  evolutions  of  the  lift 
and  drag  coefficients  imply  that  a  steady  state  is  indeed  approached. 

Figure  3  shows  plots  of  the  time  histories  of  lift,  drag  and  lift/drag  coefficients  at  various 
values  of  speed  ratios  (0  <  a  <  3.25)  and  for  time  in  the  interval  0  <  t  <  24.  As  seen  in  Fig- 
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ure  3(a),  when  the  speed  ratio  is  increased  to  2.07,  the  lift  increases  timewise  proportionally. 
However,  as  the  speed  ratio  further  increases,  lift  appears  to  initially  decrease  then  ineresises 
gradually  at  later  times.  Not  surprisingly,  the  maximum  value  of  Cl  that  can  be  achieved  by 
rotation  is  also  higher  as  the  speed  ratio  grows.  It  is  also  observed  that,  at  speed  ratios  lower 
than  2,  the  respective  lift  curves  exhibit  a  well  established  periodic  evolution.  However,  in 
the  range  of  a  >  2,  it  is  not  known  whether  the  nature  of  this  periodicity  will  continue 
if  the  time  of  investigation  is  expanded.  Apparently,  as  can  be  seen  from  Figure  3(a),  the 
cylinder  rotation  (as  a  boundary  moving  mechanism)  does  yield  a  substantial  lift  enhance¬ 
ment  as  indicated  by  the  experimental  work  of  Modi’s  group  for  airfoil/rotating  cylinders 
configurations. 

As  illustrated  by  the  drag  curve  in  Figure  3(b),  there  is  a  substantial  increase  in  drag 
when  the  speed  ratio  is  increased.  In  all  cases  considered  here,  these  drag  curves  seem  to 
converge  after  a  certain  time  and  then  oscillate  under  different  amplitudes  and  frequencies 
thereafter.  Detailed  numerical  results  on  the  effect  of  the  speed  ratio  to  the  resulting  lift/drag 
curve  are  shown  in  Figure  3(c).  In  the  range  0  <  a  <  2.07,  the  lift/drag  performance  appears 
to  improve  timewise  (for  0  <  t  <  24)  with  an  increase  of  a.  If  a  comparison  is  made  between 
a  =  2.07  and  a  =  0.05,  a  noticeable  improvement  of  the  lift/drag  performance  is  observed. 
Although  a  higher  lift/drag  ratio  is  achieved  by  increasing  the  rotation  rate  in  this  range,  the 
question  arises  whether  any  further  increase  of  a  will  result  in  a  continued  improvement  of 
the  lift/drag  ratio.  Intuitively,  it  is  natural  to  expect  a  monotonica!  increase  in  the  lift/drag 
ratio  as  a  increases  to  a  =  3.25.  However,  this  is  not  the  case  as  a  comparison  is  made 
between  a  =  3.25  and  a  —  2.07.  In  fact,  the  lift/drag  curves  illustrate  a  gradually  decrease 
in  performance  over  certain  time  interval  when  the  speed  ratio  increases  beyond  2.  Moreover, 
this  tendency  toward  lower  lift/drag  ratio  becomes  noticeable  when  a  reaches  the  highest 
value  (a  =  3.25)  considered  here.  Nevertheless,  for  all  a  considered  here,  a  significant 
increase  in  the  maximum  value  of  Cl/Cd  can  be  obtained  by  increasing  a.  However,  if  'c 
found  that  it  will  reach  its  maximum  value  at  a  much  later  time  for  higher  values  of  a. 
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3.2  Time-Periodic  Rotation  vs.  Time-Harmonic  Rotary  Oscillation 

The  previous  results  only  applied  to  constant  rotation  rates.  In  this  section  we  consider 
time-varying  rotations.  Although  the  maximum  value  of  lift/drag  ratio  can  be  achieved  at 
a  certain  optimal  constant  rotation  rate  described  above,  this  does  not  immediately  imply 
that  the  constant  rotation  will  always  yield  a  global  maximum  value  of  the  lift/drag  ratio 
over  a  specific  time  interval.  In  fact,  the  general  problem  is  to  determine  the  optimal  control 
input  among  all  time  dependent  rotation  rates.  As  a  first  step  towards  solving  this  problem, 
one  needs  to  extend  the  computational  calculations  to  the  case  of  arbitrary  time-dependent 
rotation  rates.  Since  we  shall  not  attempt  to  solve  the  necessary  conditions  corresponding  to 
the  optimal  rotation  rate  problem  and  because  the  main  goal  of  the  paper  is  to  gain  insight 
into  the  possible  form  of  an  optimal  controller,  we  shall  concentrate  on  two  periodic  inputs 
and  investigate  the  impact  of  each  control  on  this  problem.  Hence,  to  keep  the  paper  of 
reasonable  length  and  the  simulations  simple  we  restrict  our  study  to  two  periodic  inputs. 

It  is  well  known  that  when  a  cylinder  oscillates  in  a  uniform  flow,  the  associated  forcing 
oscillating  frequency  and  amplitude  can  influence  the  vortex  formulation  and  forces  response 
substantially  [42,  44].  It  has  been  experimentally  shown  that  at  Re  =  200,  the  natural 
Strouhal  frequency  of  a  non-rotating  circular  cylinder  (a  =  0)  is  approximately  Fn  —  0.185 
[43].  It  is  of  important  to  study  the  behavior  of  fluctuating  forces  at  imposed  forcing  fre¬ 
quencies  which  lie  in  the  neighborhood  of  the  natural  frequency.  The  temporal  evolutions  of 
lift,  drag  and  lift/drag  are  shown  separately  in  Figures  4(a,b,c)  for  a  time-periodic  rotation 
a(t)  =  |sin0.25<|  and  a  time-harmonic  rotary  oscillation  a(t)  =  sin0.5t,  respectively.  In  the 
case  of  time-periodic  rotation,  the  cylinder  under  control  is  rotated  in  the  counterclockwise 
direction  about  its  axis  with  a  time-periodic  angular  velocity.  Notice  that  these  two  types 
of  rotation  are  employed  by  the  same  forcing  frequency  (i.e.  F  =  0.16)  which  lies  in  the 
neighborhood  of  the  natural  frequency.  The  numerical  results  clearly  confirm  the  expected 
benefit  of  this  time-periodic  rotation  for  both  lift  and  drag  forces,  as  shown  in  Figure  4. 

In  comparing  these  two  types  of  rotation,  it  should  be  noted  that  rotating  in  the  same 
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direction  causes  the  lift  curve  * "  be  shifted  upwards  due  to  the  nature  of  rotation,  while  the 
drag  curve  is  shifted  downwards.  In  terms  of  performance,  this  corresponds  to  an  increase  of 
the  time-averaged  lift  force  in  the  time-span  of  the  investigation,  while  in  the  same  time  in¬ 
terval,  a  substantial  reduction  of  the  time-averaged  drag  as  well.  The  resulting  improvement 
of  the  lift/drag  ratio  is  shown  in  Figure  4(c). 

To  demonstrate  the  influence  of  this  time-periodic  rotation  on  the  temporal  development 
of  these  force  coefficients,  two  additional  values  of  the  forcing  frequency  were  tested.  In 
Figure  5,  all  the  time  histories  of  lift,  drag  and  lift/drag  coefficients  are  shown  for  a  forcing 
frequency  (i.e.  F  —  0.08)  which  is  1/2  the  frequency  in  Figure  4.  As  it  can  be  seen  from 
Figure  5(a),  except  in  an  initial  stage,  lift  is  increased  for  the  time-periodic  rotation  (a(t)  = 
|  sin  0.125<|)  when  compared  to  the  time-harmonic  rotary  oscillation  (o(t)  =  sin0.25f).  Also 
the  drag  is  reduced  as  illustrated  in  Figure  5(b).  However,  the  amount  of  drag  reduction  is 
not  significant  over  the  time-span  of  the  investigation.  As  would  be  expected,  the  resulting 
lift/drag  curves  shown  in  Figure  5(c)  exhibit  almost  the  same  behavior  as  those  lift  curves 
illustrated  in  Figure  5(a). 

Similar  results  are  shown  in  Figure  6,  where  a  higher  rotation  frequency  (i.e.  F  =  0.32)  is 
imparted  to  the  cylinder.  As  before,  the  time-periodic  rotation  at  this  frequency  enhances  the 
lift  performance  when  compared  to  the  time-harmonic  rotary  oscillation.  However,  notice 
that  the  corresponding  drag  curves  are  oscillating  about  an  “average  value”  in  the  time 
interval  under  consideration,  for  both  types  of  rotation.  The  improvement  in  the  lift/drag 
curve  shown  in  Figure  6(c)  is  more  noticeable  when  compared  to  the  lift/drag  curve  in 
Figure  4(c).  Notice  that  none  of  the  frequencies  considered  in  Figures  5  and  6  are  in  the 
neighborhood  of  the  natural  frequency. 

In  a  detailed  examination  of  all  comparisons  between  two  types  of  rotations  mentioned 
above,  it  is  found  that  lift  enhancement  is  not  always  timewise  in  the  time  interval  under 
investigation.  The  exact  time  at  which  the  time-periodic  rotation  outperform  the  time- 
harmonic  rotary  oscillation  depends  on  the  forcing  frequency  imparted  to  the  cylinder.  This 
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occurs  earlier  for  a  higher  value  of  forcing  frequency.  However,  the  drag  curves  have  a  signif¬ 
icantly  different  character.  The  decrease  or  increase  in  drag  depends  on  whether  its  forcing 
frequency  lies  in  the  neighborhood  of  natural  frequency  or  not.  If  the  forcing  frequency  is  in 
the  neighborhood  of  the  natural  frequency  as  the  case  in  Figure  4(b),  the  drag  curve  exhibits 
a  substantial  reduction.  On  the  other  hand,  for  the  forcing  frequencies  which  are  not  in  the 
neighborhood  of  the  natural  frequency  as  shown  in  Figures  5(b)  and  6(b),  the  drag  curves 
changes  only  slightly.  In  addition,  for  these  forcing  frequencies,  the  respective  lift/drag 
curves  are  similar  to  the  corresponding  lift  curves.  This  is  due  to  the  minor  variations  in  the 
drag  curves. 

3.3  Effects  of  Forcing  Frequency  and  Angular  Amplitude 

Because  of  the  above  observations,  it  is  worthwhile  to  investigate  the  effects  of  vari¬ 
ous  control  parameters  upon  the  force  coefficients.  A  comparisons  of  force  coefficients  are 
shown  in  Figure  7,  corresponding  to  all  three  forcing  frequencies  described  above.  It  can  be 
readily  seen  that  these  forcing  frequencies  have  considerable  influence  on  the  amplitude  and 
frequency  of  the  oscillatory  forces. 

In  terms  of  performance,  Figure  7(a)  presents  the  evolution  of  lift  coefficients  and  shows 
no  clear  advantage  of  changing  the  forcing  frequency.  However,  as  shown  in  Figures  7(a,b), 
the  rotation  ( a(t )  =  |sin0.25t|)  which  lies  in  the  neighborhood  of  the  nature  frequency 
achieves  a  higher  value  of  ( Ci)max ,  and  yields  a  slightly  larger  value  of  drag  when  compared 
to  the  other  two  cases  (namely,  a(t)  =  |sin0.125t|  and  <*(f)  =  |sin0.5f|).  The  lift/drag 
curves  shown  in  Figure  7(c)  exhibit  similar  temporal  evolutions. 

It  is  also  interesting  to  study  the  effect  of  angular  amplitude  on  the  temporal  evolution 
of  forces  while  the  forcing  frequency  is  fixed  as  a  constant.  Figure  8  shows  that  resulting 
forces  on  the  cylinder  can  differ  significantly  at  different  angular  amplitudes  for  a(t)  = 
A\  sin0.314£|.  This  rotation  corresponds  to  a  forcing  Strouhal  number  of  0.2  which  is  in  the 
neighborhood  of  the  natural  Strouhal  number  of  0.185.  The  angular  amplitudes  considered 
were  A  =  1.0,2.07  and  3.25.  Apparently,  as  can  be  seen  from  these  figures,  a  larger  angular 
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amplitude  definitely  yields  an  incremental  lift  coefficient  over  the  time-span  of  investigation 
(0  <  t  <  36).  However,  initially  the  drag  increases  with  an  increase  of  A ,  then  after  a  certain 
time  it  oscillates  with  almost  the  same  amplitude  and  frequency  around  an  averaged  value. 
Consequently,  this  leads  to  a  substantial  improvement  in  lift/drag  with  increasing  .4,  as 
clearly  shown  in  Figure  8(c).  Nevertheless,  the  effect  of  angular  amplitude  is  very  noticeable 
when  compared  to  the  effect  of  the  forcing  frequency  shown  in  Figure  7. 

3.4  Time- Averaged  Value  of  Forces 

We  also  examined  the  averaged  values  of  the  force  coefficients  over  the  time  span  of 
investigation  as  various  control  parameters  are  altered.  The  time-harmonic  rotary  oscillation 
(a(t)  =  sin7rF<)  described  in  §3.2  was  considered  with  three  different  values  of  forcing 
frequency.  The  forces  were  averaged  with  respect  to  the  time  (0  <  t  <  24).  Figure  9(a) 
shows  the  time-averaged  values  of  lift,  drag,  lift/drag  as  the  forcing  frequency  F  is  varied 
between  0.08  and  0.32.  This  figure  shows  that  the  forcing  frequency  has  a  considerable 
influence  on  the  time-averaged  values  of  the  force  coefficients.  The  local  maximum  values  of 
tiine-averaged  lift,  drag  and  lift/drag  ratios  correspond  to  the  forcing  frequency  which  lies 
in  the  neighborhood  of  the  natural  frequency.  This  particular  feature  was  also  observed  in 
the  numerical  results  of  Mo  [23]  where  it  was  shown  that  the  drag  peak  occurs  at  the  forcing 
frequency  equal  to  the  natural  frequency. 

As  for  the  cases  of  time-periodic  rotation,  variations  of  time-averaged  forces  coefficients 
with  respect  to  the  forcing  frequency  are  presented  in  Figure  9(b).  As  illustrated  in  the 
figure,  a  forcing  frequency  in  this  range  (i.e.  0.08  <  F  <  0.32)  has  little  effect  on  the 
time-averaged  forces.  Although  the  difference  in  time-averaged  drag  is  minor,  the  forcing 
frequency  which  lies  in  the  neighborhood  of  the  natural  frequency  (F  =  0.185)  corresponds 
to  a  larger  time-averaged  drag  and  a  smaller  time-averaged  lift. 

The  effect  of  angular  amplitude  on  the  time-averaged  values  of  lift,  drag  and  lift/drag 
coefficients  is  shown  in  Figure  9(c),  for  a(t)  =  A|sin0.314f|  averaged  over  0  <  t  <  36.  For 
A  in  the  range,  1  <  A  <  3.25,  all  the  time-averaged  values  are  almost  linearly  proportional 
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to  the  angular  amplitude.  Significant  increases  in  lift  coefficients  with  increasing  angular 
amplitude  is  particularly  noticeable.  This  can  be  demonstrated  by  comparing  the  case  A  — 
3.25  with  .4=1.  It  represents  a  240%  increment  of  lift  performance.  However,  a  slight 
increment  in  the  drag  coefficients  with  increasing  angular  amplitude  is  also  observed.  A 
moderate  improvement  of  time-averaged  lift/drag  ratio  is  also  seen. 

In  the  case  of  a  constant  speed  of  rotation,  it  is  also  worthwhile  to  study  the  effect  of 
speed  ratios  on  the  time-averaged  lift,  drag  and  lift/drag  coefficients.  These  results  are 
shown  in  Figure  10(a)  for  a  in  the  range,  0  <  a  <  3.25,  and  in  the  time  interval  0  <  t  <  24. 
It  illustrates  that  the  time-averaged  lift  is  almost  linearly  proportional  to  the  speed  ratio, 
while  the  time-averaged  drag  remains  as  a  constant  value  up  to  «  =  2,  then  monotonically 
increases  with  speed  ratio  thereafter.  Most  importantly,  the  time-averaged  lift/drag  is  not 
linearly  proportional  to  the  speed  ratio.  As  shown  in  the  figure,  the  highest  value  of  the 
speed  ratio  a  =  3.25  considered  here  is  not  the  optimal  constant  rotation  rate  corresponding 
to  the  maximum  value  of  time-averaged  lift/drag.  The  maximum  value  occurs  at  a  lower 
speed  ratio,  approximately  a  =  2.38,  and  it  represents  a  substantial  increase  of  440%  over 
the  lower  speed  ratio  a  =  0.5.  In  Figure  10(b),  the  variation  of  the  (total  lift) /(total  drag) 
force  ratio  with  respect  to  the  speed  ratio  is  shown  for  a  in  the  range  0  <  a  <  3.25.  Although 
a  maximum  is  achieved  at  a  value  between  a  =  2.0  and  a  =  2.38,  it  should  be  noted  that 
this  optimal  speed  ratio  is  not  necessarily  the  same  optimal  value  as  described  in  Figure 
10(a). 

The  results  presented  in  Figures  10(a,b)  demonstrate  an  effective  way  of  improving 
lift/drag  performance  by  changing  the  rotation  rate  and  illustrate  the  important  of  selecting 
a  proper  rotation  rate  in  order  to  maximize  performance.  If  one  formulate  this  problem  as 
an  optimal  control  problem,  then  the  time-averaged  lift/drag  ratio  in  Figure  10(a)  represents 
the  cost  function  and  the  goal  is  to  find  an  optimal  control  a,,  among  a  set  of  restricted  con¬ 
trol  parameters  (constant  values),  that  will  provide  the  maximum  value  of  the  time-averaged 
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performance  functional 
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where  Tf  is  the  final  time.  Similarly,  the  curve  in  Figure  10(b)  represents  the  optimal  control 
problem  defined  by  maxim  '  dng  the  cost  functional 
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In  theory,  these  optimal  control  problems  can  be  solved  by  applying  necessary  conditions 
for  distributed  parameter  systems.  However,  computational  methods  for  such  necessary 
conditions  are  complex  and  not  yet  fully  developed.  The  results  presented  here  may  be  used 
to  guide  and  test  future  computational  schemes  based  on  optimal  control  theory. 


3.5  Synchronizations  of  Cylinder  and  Wake 


The  synchronization  of  cylinder  and  wake  has  long  been  known  to  be  an  important 
component  of  vortex-induced  oscillations  [33].  A  detailed  study  of  various  types  synchro¬ 
nization  for  a  body  oscillating  transversely  in  a  uniform  stream  can  be  found  in  Williamson 
and  Roshko  [44],  For  the  case  of  time-periodic  rotatior.  considered  here,  it  is  natural  to  ask 
whether  such  synchronization  can  occur  and  how  well  the  numerical  results  can  predict  the 
occurrence  of  this  important  phenomenon.  To  the  best  of  our  knowledge,  the  current  study 
is  the  first  work  to  investigate  synchronization  under  the  particular  form  of  time-periodic 
rotation  described  in  previous  sections. 

An  examination  of  the  responses  in  Figure  4,  shows  that  the  combined  system  of  cylinder 
and  wake  will  be  “locked  in”  by  an  imposed  forcing  frequency.  This  synchronization  of 
the  cylinder  and  wake  is  due  to  the  fact  that  the  forcing  frequency  of  rotation  ( F  =  0.16) 
lies  in  the  neighborhood  of  the  natural  frequency  (F„  =  0.185).  Notice  that  in  the  case 
of  time-periodic  rotation  shown  in  Figure  4,  both  lift  and  drag  curves  oscillate  with  the 
forcing  frequency  (corresponding  to  a  time  period  of  T  —  12.5),  clearly  exhibiting  a  periodic 
response.  However,  in  the  case  of  time-harmonic  rotary  oscillation,  the  lift  curve  oscillates 
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with  the  same  forcing  frequency  ( T  =  12.5)  while  drag  curve  oscillates  with  the  period  of 
Tf  2.  Consequently,  the  lift/drag  ratios  oscillate  at  the  same  frequency  ( T  =  12.5)  for  both 
types  of  rotation. 

For  the  case  of  time-periodic  rotation  a(t )  =  A\  sin0.314£|,  we  extend  our  observation  to 
a  relatively  longer  time.  For  0  <  t  <  36,  an  examination  of  these  force  curves  for  A  —  1.0 
in  Figure  8  exhibits  a  periodic  response  with  a  frequency  ( F  =  0.2)  precisely  equal  to  the 
input  forcing  frequency  (i.e.  T  =  10).  Although  this  periodic  behavior  is  not  established  for 
A  =  2.07  and  3.25,  the  corresponding  curves  are  almost  periodic  in  time.  In  order  to  confirm 
this  periodicity,  a  sequence  of  instantaneous  streamlines  are  shown  in  Figure  11.  In  Figure 
11,  each  plot  is  separated  with  an  interval  of  one  time  period.  These  streamlines  are  plotted 
in  a  frame  fixed  with  the  undisturbed  fluid.  The  periodicity  of  the  flow  is  clearly  noticeable. 
Two  opposite-sign  vortices  are  shed  alternately  on  opposite  sides  of  the  cylinder  at  each  cycle 
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of  rotation.  The  vortex  formation  in  the  wake  is  similar  to  the  case  of  a  non-rotating  cylinder 
(a  =  0).  However,  the  midline  of  the  vortex  street  has  been  displaced  slightly  upwards  due 
to  the  nature  of  rotation  (in  the  counterclockwise  direction). 

In  order  to  identify  the  range  of  frequency  for  this  fundamental  synchronization,  we  im¬ 
pose  a  rotation  rate  a(t)  =  |sin0.283f|  with  a  forcing  frequency  ( F  =  0.18)  which  is  the 
neighborhood  of  the  natural  frequency.  The  time  histories  of  lift,  drag  and  lift/drag  coeffi¬ 
cients  shown  in  Figures  12(a,b,c)  clearly  demonstrate  the  periodic  behavior  of  the  response. 
The  corresponding  streamline  plots  (at  each  instant)  are  presented  in  Figure  13.  These  results 
show  that  there  exists  a  range  of  forcing  frequencies  in  which  fundamental  synchronization 
will  occur.  However,  the  precise  range  of  forcing  frequencies  leading  to  synchronization  has 
not  yet  been  determined  (computational  time  is  the  limiting  factor). 

In  the  case  of  time-harmonic  rotary  oscillation,  the  effects  of  the  forcing  frequency  and 
amplitude  on  a  cylinder  wake  have  been  investigated  experimentally  by  Tokumaru  and  Di- 
motakis  [40].  Several  vortex  formations  were  observed  in  the  wake.  Their  experiments  dealt 
with  a  range  of  amplitudes  and  frequencies  at  a  Reynolds  number  of  Re  —  1.5  x  104,  By  fixing 
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the  reduced  amplitude  A  in  their  experiments,  four  qualitatively  different  vortex  shedding 
modes  were  identified  when  the  forcing  frequency  was  increased.  Although  their  experiments 
were  conducted  at  a  Reynolds  number  higher  than  the  current  study  (Re  =  200).  we  con¬ 
ducted  a  similar  investigation.  Here,  a  particular  value  of  the  rotation  rate  a (t)  —  2. sin  3.1-}/ 
(corresponding  to  A  =  2  and  F  =  1.0)  was  tested.  Figure  14  shows  the  lift  and  drag  histories 
up  to  t  =  36.  These  curves  are  clearly  periodic  in  nature.  Under  this  forcing  frequency,  t  he 
lift  curve  oscillates  with  a  period  of  T  =  10,  while  the  drag  curve  oscillates  with  a  period 
of  T  =  5.  The  instantaneous  streamlines  plots  are  presented  in  Figure  15,  and  they  show 
a  time  periodic  flow  pattern.  Moreover,  these  results  indicate  that  rotation  may  provide  an 
effective  control  of  the  cylinder  wake. 

4.  CONCLUSIONS 

An  algorithm  for  computing  the  viscous  flow  past  a  rotating  cylinder  is  presented  and 
applied  to  the  problem  of  controlling  cylinder  forces  by  rotation.  Several  fundamental  types  of 
rotation  were  considered.  Using  time-periodic  rotations  leads  to  a  considerable  improvement 
in  the  force  coefficients  and  was  shown  to  be  very  effective,  especially  compared  to  time- 
harmonic  rotary  oscillations.  These  results  are  significant  because  they  show  a  proper  choice 
of  the  rotation  rate  can  lead  to  improved  flow  fields.  Very  precise  periodicity  of  the  force 
for  certain  cases  was  established,  and  this  periodic  behavior  has  considerable  impact,  on 
controlling  the  vortex  formation  in  the  cylinder  wrake.  For  the  case  of  a  constant  speed  of 
rotation,  two  optimal  control  problems  were  considered  and  solved  computationally. 

These  results  demonstrate  the  feasibility  of  using  boundary  mechanisms  for  controlling 
unsteady  flows,  and  consequently  can  be  applied  to  enhance  the  performance.  Using  such 
mechanisms  as  a  controller  allows  us  to  formulate  a  wide  variety  of  optimal  control  problems 
for  fluid  flow  systems.  Modifications  of  existing  numerical  algorithms  needed  for  such  cont  rol 
problems  depend  on  performance  and  design  constraints.  For  example,  one  may  need  the 
maximum  (or  minimum)  lift  to  drag  ratio  in  order  to  sustain  a  particular  maneuver  of  a 
supermaneuverable  aircraft.  Because  of  the  complexity  and  importance  of  the  relationship 
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between  vortex  motion  and  cylinder  forces,  the  first  step  in  control  design  may  be  to  seek  a 
specific  type  of  rotation  rate  that  will  match  all  proposed  goals. 

A  precise  understanding  of  time-dependent  moving  surfaces  in  boundary  layer  control 
may  provide  an  effective  way  for  lift  enhancement  and  drag  reduction.  By  treating  the 
rotation  rate  as  a  control  variable  in  this  model,  we  will  eventually  be  interested  in  finding 
the  optimal  control  (i.e.  a  time  history  of  the  rotation  rate)  that  maximizes  (or  minimum) 
the  lift-to-drag  ratio  over  a  fixed  time  interval.  Although  here  the  optimal  control  problem 
associated  with  the  constant  rotation  rate  was  solved  by  direct  computations,  it  is  still 
important  to  explore  the  possible  implementation  of  a  computational  algorithm  to  calculate 
the  optimal  solution  for  the  more  general  problems.  The  tools  developed  here  can  be  used 
to  investigate  fundamental  questions  regarding  control  of  separated  flows  by  using  various 
boundary  control  mechanisms.  Future  work  needs  to  be  done  in  the  development  of  new 
computational  algorithms  for  solving  complex  optimal  flow  control  problems. 
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Re  =  200  with  various  constant  speed  ratios  (0.05  <  a  <  3.25). 
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FIGURE  5:  Comparison  of  temporal  evolution  of  the  lift  (a),  drag  (b)  and  lift/drag  (c) 
coefficients  for  a  time-periodic  rotation  a(t)  =  |sin0.125£|  with  a  time-harmonic  rotary 
oscillation  a(t)  =  sin0.25f  at  Re  =  200  for  0  <  t  <  24. 


FIGURE  6:  Comparison  of  temporal  evolution  of  the  lift  (a),  drag  (b)  and  lift/drag  (c)  co 
efficients  for  a  time-periodic  rotation  a(t)  =  |  sin  0.5t|  with  a  time-harmonic  rotary  oscillatioi 
a(t )  =  sint  at  Re  =  200  for  0  <  t  <  24. 


FIGURE  8:  Temporal  evolution  of  the  lift  (a),  drag  (b)  and  lift/drag  (c)  coefficients 
for  a  time-periodic  rotation  a(t)  =  A|  sin  0.314f)  at  Re  -  200  with  various  amplitudes  of 
A  =  1.0,2.07  and  3.25  for  0  <  t  <  36. 
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FIGURE  9:  Variation  of  time-averaged  forces  coefficients  with  respect  to  the  forcing 
frequency  and  the  angular  amplitude:  (a)  a(t)  =  sin  nFt  and  0.08  <  F  <  0.32;  (b) 
a(t)  =  |sinir(F/2)t|  and  0.08  <  F  <  0.32;  (c)  a(t)  =  A\  sin0.314f|  and  1  <  A  <  3.25. 
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FIGURE  10:  Effect  of  the  speed  ratio  on  time-averaged  lift,  drag  and  lift/drag  coefficients 
(a)  and  on  total  lift/total  drag  force  ratio  (b)  for  0  <  a  <  3.25. 


40 


FIGURE  11:  Instantaneous  streamlines  plots  for  Re  =  200,  a{t)  —  | sin 0.31 4fj  (F  =  0.2), 
viewed  from  a  frame  fixed  with  the  undisturbed  fluid,  (a)  t  =  16,  (b)  t  =  26,  (c)  t  =  36. 
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FIGURE  13:  Instantaneo 
viewed  from  a  frame  fixed  \ 


FIGURE  15:  Instantaneous  streamlines  plots  for  Re  =  200,  a(t )  =  2sin3.14f  ( F  =  l.( 
and  A  —  2),  viewed  from  a  frame  fixed  with  the  undisturbed  fluid,  (a)  t  =  1G,  (b)  t  =  26 
(<:}  t  =  3G. 
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